Flow Between Two Sites on a Percolation Cluster 
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We study the flow of fluid in porous media in dimensions d = 2 and 3. The medium is modeled by 
bond percolation on a lattice of L'' sites, while the flow front is modeled by tracer particles driven 
by a pressure difference between two fixed sites ("wells") separated by Euclidean distance r. We 
investigate the distribution function of the shortest path connecting the two sites, and propose a 
scaling Ansatz that accounts for the dependence of this distribution (i) on the size of the system, 
L, and (ii) on the bond occupancy probability, p. We confirm by extensive simulations that the 
Ansatz holds for d = 2 and 3, and calculate the relevant scaling parameters. We also study two 
dynamical quantities: the minimal traveling time of a tracer particle between the wells and the 
length of the path corresponding to the minimal traveling time "fastest path" , which is not identical 
to the shortest path. A scaling Ansatz for these dynamical quantities also includes the effect of 
finite system size L and off-critical bond occupation probability p. We find that the scaling form 
for the distribution functions for these dynamical quantities for d = 2 and 3 is similar to that 
for the shortest path but with different critical exponents. The scaling form is represented as the 
product of a power law and three exponential cutoff functions. We summarize our results in a table 
which contains estimates for all parameters which characterize the scaling form for the shortest path 
and the minimal traveling time in 2 and 3 dimensions; these parameters are the fractal dimension, 
the power law exponent, and the constants and exponents that characterize the exponential cutoff 
functions. 



I. INTRODUCTION 



Percolation theory is a paradigmatic model for connec- 
tivity, originally introduced as a mathematical subject in 
the late 1950s [|l|. Thereafter, percolation theory has 
been found useful to characterize many disordered sys- 
tems [|-|. 

The simplest percolation model is a lattice of bonds 
occupied with probability p. Neighboring bonds are con- 
sidered to be connected if both are occupied. A set of 
sites connected by bonds is called a cluster. As p in- 
creases, new clusters are formed and previously existing 
clusters not only grow, but become connected as more 
sites are occupied in the system. At a critical value of 
p, Pc (known as the percolation threshold), one spanning 
cluster appears and provides overall connectivity. Just at 
the critical point, the incipient infinite percolation cluster 
is an example of a random fractal that is a useful model 
for real disordered systems. While the actual thresh- 
old value Pc depends on the particular lattice chosen, 
the behavior of the properties measured near the perco- 
lation threshold is universal, depending only on the di- 
mensionality of the system. For example, the mass M of 
a cluster diverges at the percolation threshold as a power 
M - |p-pcr^ where 7 = 43/18 {d = 2) and 1.795±0.005 



{d — 3) regardless of the type of lattice [Q-^. This uni- 
versality concept is extremely useful as it means that we 
can understand a new system knowing that its critical 
exponents are the same as previously-studied systems. 
Also, off-lattice continuum percolation appears to have 
the same exponent values as those for lattice percolation 
i|. 

A comprehensive set of exactly- and numerically- 
calculated critical exponents is now available to describe 
many of the features of percolation and the perco- 
lation paradigm has been applied to problems of practi- 
cal interest in heterogeneous chemistry Q], polymer sci- 
ence pO| and transport phenomena in disordered systems 
pT| . In studying transport phenomena in disordered sys- 
tems, one must couple additional processes to the geo- 
metrical features of the percolation representation. Typ- 
ical examples include the problems of diffusion, reaction 
and flow through porous media, as well as the metal- 
to-insulator transition in polymer composites and alloys. 
These are systems in which the interplay between struc- 
ture and phenomenology must be investigated in detail, 
as it might be the relevant factor determining optimum 
material properties. 

A useful concept — first put forward by Ambegaokar, 
Halperin and Langer |l^ — is that electrical transport in 
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disordered media with a broad distribution of conduc- 
tance values is dominated by those regions where the 
conductances are larger than some critical value a^. This 
critical value is the largest conductance such that the set 
of conductances above this threshold value still preserves 
the global connectivity of the system. In percolation ter- 
minology, this cluster would be analogous to the conduct- 
ing spanning cluster. It is then possible to reduce the 
general problem of transport in a highly connected dis- 
ordered media with a broad distribution of conductances 
to a percolation problem at criticality. Once the "criti- 
cal cluster" is identified within the disordered geometry, 
one can estimate macroscopic transport properties of the 
system. This approximation, commonly called "the crit- 
ical path method" has been extended by Katz and 
Thompson jist to estimate transport properties (e.g., 
permeability and electrical conductivity) in disordered 
materials. Thus, although the majority of studies of fluid 
flow through disordered media are for systems well above 
the threshold (e.g., sandstone), if the distribution of per- 
meabilities is sufficiently broad, we can still make use of 
percolation concepts to model the relevant geometry of 
well-connected disordered media. Therefore, percolation 
theory is certainly an appropriate description for a large 
number of disordered systems (see also Ref. and ref- 
erences therein). This is the basis of our study and also 
the basis of the vast literature on fluid flow in percolation 
clusters. 

The aim of the present paper is to discuss the potential 
application of percolation theory as a convenient geomet- 
rical model for understanding numerous aspects of flow 
through porous rocks |l^,|l^. Special emphasis will be 
given to the study of oil displacement, i.e., how hydrocar- 
bons propagate through geological formations between a 
pair of wells in the oil field. This work could also be 
applied to the breakthrough time for contamination of 
a water supply, or the time for released radioactive ma- 
terial to get from a leaking nuclear repository into the 
biosphere. 

Oil fields are extremely complex, containing geological 
heterogeneities on a wide range of length scales from cen- 
timeters to kilometers. These heterogeneities, caused by 
the sedimentary processes that deposited the rocks and 
the subsequent actions on the rock, such as fracturing by 
tectonic forces and mineral deposition from aquifer flow, 
have a significant impact on hydrocarbon recovery. 

The most common method of oil recovery is by dis- 
placement. Either water or a miscible gas (carbon diox- 
ide or methane) is injected in a well (or wells) to displace 
the oil to other wells. Ultimately the injected fluid will 
break through into a production well where it must be 
separated from the oil, which is a very costly process. 
Once the injected fluid has broken through, the rate of 
oil production declines as more injected fluid is produced. 
For economic purposes it is important to know when the 
injected fluid will break through and what the rate of de- 
cline of oil production will be so that the economic limit 
of production can be determined. 



Because the sedimentary process that produces the 
porous rock through which the fluid flows is very chaotic, 
the rock is highly heterogeneous. However, in many cases 
the rock can be separated into two types — high perme- 
ability ( "good" ) and low or zero permeability ( "bad" ) — 
and for all practical purposes we can assume the flow 
takes place only in the good rock. The spatial distribu- 
tion of the rock types is often close to random, in which 
case the classical percolation problem is a good approx- 
imation. The place of the occupancy probability p is 
taken by the volume fraction of the good rock, called the 
net-to-gross ratio in the petroleum literature. 

We have very little direct knowledge about the spatial 
distribution of rock properties in a reservoir. Direct mea- 
surements arc limited to samples that represent a fraction 
of « 10^ ^■^ of the total reservoir volume. These samples 
are taken from the well locations. Elsewhere the prop- 
erties have to be inferred from knowledge of the general 
geological environment and by analogy with other reser- 
voirs or surface outcrops. Hence, there is a great deal of 
uncertainty in our prediction of the spatial distribution 
of these rock properties. This leads to an uncertainty in 
our ability to predict the flow performance, principally 
the breakthrough time and the production decline rate. 
We need to estimate the uncertainty accurately so that 
economic risk evaluations can be made. 

The conventional approach to this problem is to build a 
detailed (numerical) model of the rock properties. These 
models will honor the one and two point statistics ob- 
served from the wells and analogue outcrops. The models 
must also agree with the observed data values at wells. 
These models are statistical in nature and conventionally 
one samples realizations from the models and performs 
numerical flow calculations on the realizations to give a 
Monte Carlo prediction of breakthrough and production 
decline. Unfortunately this process is so time consuming 
as to be impractical in many cases. Typically the flow 
simulation can take several hours on reasonable worksta- 
tions. When hundreds of realizations are sampled to get 
good statistics, the total computing time becomes un- 
wieldy. Thus there is a strong need to make this more 
efficient so as to come up with very quick, but accurate, 
predictions of recovery and the uncertainty due to the 
lack of knowledge of the underlying rock properties. 

The purpose of this work is to accomplish the goals 
described above using methods derived from percolation 
theory. It is based on two key assumptions. One is that 
for many cases the permeability disorder can be approx- 
imated by either permeable or impermeable rock. For 
example the reservoir may have been deposited by me- 
andering river belts in which case the good sand occurs 
as "packages" in an insulating background. The other 
assumption is that the fiow paths are strongly controlled 
by the permeability disorder and not strongly modified 
by the fiow dynamics themselves. Again there are many 
cases when this is reasonable, in particular if the viscos- 
ity ratio between the injected and displaced fiuids is not 
too large or when the system is highly disordered. 
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Under these assumptions we can consider the under- 
lying heterogeneity to be that of a percolating system 
(not necessarily at threshold). This has previously been 
done to study the fraction of sand connected to well pairs 
prf . We then look at the dynamic displacement on this 
cluster where the flow is controlled by Darcy's law (analo- 
gous to Ohm's law in electrical current flow). We assume 
that the injected fluid can be treated as a passive tracer 
(i.e., one that is not absorbed by the rock) that is con- 
vected along these flow paths. Then the breakthrough 
time is the same as the first passage time and is strongly 
controlled by the shortest path length between the wells; 
the post breakthrough production decline is controlled by 
the longer paths. Figure |l|a illustrates a typical percola- 
tion backbone and the shortest path between two points. 
Note that the lines here represent not microscopic pores 
but rather sand bodies whose size are of the order of tens 
of meters. 

The rest of the paper is organized as follows. In Sec. II, 
because of its importance to flow in the cluster, we study 
the distribution of shortest paths between two sites in a 
percolation cluster. In Sec. Ill we study the distribution 
functions for the dynamic quantities: minimal traveling 
time, and length of the path corresponding to minimal 
traveling time. Finally, in Sec. IV, we draw some conclu- 
sions and discuss possible future work. 



II. SHORTEST PATH 

This section deals with the distribution of the shortest 
path between two sites on a percolation cluster. Because 
of the qualitative resemblance between the shortest path 
and the minimal traveling time of a tracer particle, the 
first step in understanding fluid transport between two 
sites in a percolation system is to characterize the geo- 
metrical properties of the shortest connecting path. For 
example, if we assume that the traveling time along a 
path is proportional to the path length (i.e., all velocities 
are equal), then we can obtain a rough estimate for the 
traveling time from purely geometrical arguments. 



A. Basic distribution functions 

The shortest path or chemical distance, i, between two 
sites on a percolation cluster is defined as the shortest 
path connecting the two sites (Fig. |l|) ||l^,|l^. The typi- 
cal value £* of the shortest path between two sites on a 
cluster scales with the geometrical distance, r, between 
these points as 



where 



1.13 ±0.02 [d = 2] 
1.374 ±0.02 [d = 3] 



(1) 



(2) 



is the fractal dimension of the shortest path pCl| , pi| . 

Consider a hypercubic lattice of sites. All infor- 
mation about the distribution of shortest paths is con- 
tained in the joint probability density function P{r,£), 
i.e., the probability that two sites on the same spanning 
cluster are separated by geometrical distance r and chem- 
ical path £. We sum over all chemical paths £ to calculate 
the probability distribution that the Euclidean distance 
between two sites is r, 



P{r) EE / P{r,£)d£. 



(3) 



Similarly, we obtain the probability distribution that two 
sites are separated by the chemical distance £ by sum- 
ming over all possible geometrical distances. 



P{£) = / P{r, £)dr . 



(4) 



Given that the shortest distance between these sites is 
£ , the conditional probability that the geometrical dis- 
tance between two sites is r is H 



P{r\£) = 



P(r,£) 
P{£) 



(5) 



For isotropic media this function has been studied exten- 
sively and P{r\£) is of the form 



Pir\£)^A,[^Y 



exp —a 



where 



and 



1 - £> rfmin - 1 ' 



V = l/dmin- 

For d = 2, Ziff recently argued that 

5,. - 1 = 25/24 M = 2]. 



(6) 



(7) 



(8) 



(9) 



Our simulations confirm the analytical form of the P{r\£) 
as well as these values of v and gr (see Fig. ||a). 

The function of interest to us is the conditional proba- 
bility for two sites to be separated by the shortest path £ , 
given that the geometrical distance between these sites 
is r 



P{£\r) 



P{r) 



(10) 



iFrom (|0|) and (|), we see that P{r\£) and P{£\r) are 
related 



P{£\r) = P{r\£) 



m 

P{r)- 



(11) 
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At the percolation threshold, it has been shown [|5| that, 
in analogy with (^), 



P{i\r) 



where 



1 ^ e 



-gt 



exp -a , 

/•"mm 



^_ {gr-l) + {2-df) ^ 



= vl{l -i') = 



and 



, _ r 91/48 [d= 2] 

- \ 2.524 ± 0.008 [d = 3] 



(12) 



(13) 



(14) 



(15) 



is the fractal dimension of the incipient infinite cluster 
g. Substituting (|) into (|l3|), we find for d = 2 



gi = 2.01 ±0.02 



[d = 2]. 



(16) 



The probability distribution of more practical inter- 
est is P'(^|r), defined in the same way as P{t\r) but for 
any two randomly-chosen points separated by geomet- 
rical distance r and on the same cluster, but not nec- 
essarily on the incipient infinite cluster p^ . As shown 
in Appendix ^ P'{l\r) has the same scaling form as in 
Eq. (O), but with gi replaced by 



gi = ge + 



d-d 



f 



(17) 



Figure ||b illustrates the difference between gi and g^. 



B. Distribution of shortest path 

The complete scaling form of P'(€|r), which accounts 
also for finite size effects and off-critical behavior, has 
been studied for d = 2 and reported in [l5|. Specifically, 
the following Ansatz has been proposed [|15|| 



P'{l\r) 



1 



7" inin \ 7* min 



'9i 



h 



(18) 



where C ~ |P ~ Pc| ''is the pair connectedness length, 
and the scaling functions have the form 



/i(x) = exp(-aa; '^), 



/2(a;)=exp(-6a;'^), 



(19) 



(20) 



and 



f:i{x) = exp(-cx). 



(21) 



The function /i accounts for the lower cut-off due to 
the constraint £ > r, while /2 and /a account for the up- 
per cut-offs due to the finite size effect and due to the 
finite correlation length respectively. Either /2 or /s be- 
comes irrelevant, depending on the magnitudes of L and 
^: for L < ^, /2 dominates the upper cut-off, otherwise /a 
dominates. We assume the independence of the finite size 
effect and the effect of the concentration of the occupied 
sites, so that Eq. (18) can be represented as a product of 



the terms which are responsible for the finite size effect 
(/2) and the effect of the concentration (/a). Simulations 
for d = 2 |l5| have been used to test this assumption. 



C. Shortest path in three dimensions 

1. Behavior at criticality 

Here we extend the study of P'(^|r) to d = 3. We nu- 
merically test the scaling conjecture (|l^) exactly at the 
percolation threshold p = Pc — in which case ^ = oo so 
/a = /(O) — const. We build clusters using the Leath 
algorithm jl^Jl^,|2^. Since the Leath algorithm corre- 
sponds to the process of selecting a random point on the 
lattice, the probability P'{£\r) is equal to the probabil- 
ity that a pair of randomly selected points has chemical 
distance £ and geometrical distance r, given that they be- 
long to the same cluster, a cluster that is not necessarily 
the infinite cluster. Hence Eq. (nsh reduces to 



P\£\r) 



1 



e 



-9 1 



fl 



(22) 



Figure l^a shows that, in the range r'^mm < £ < L'^min, 
P'{£\r) has power-law behavior with slope 



5; = 2.3 ±0.1, 



[d = 3] 



(23) 



and rapidly vanishes for £ < r'^niin and for £ > L'^min. 
To determine the functions /i and /2, we compute the 
rescaled probability distribution 



P'{£\r) £^'<^ 7.-''min(9«-l)^ 



(24) 



and plot it against scaling variable x ; 
Fig. 0b) using the value dmin = 1.374. 



$(x)=A/i(x)/2 



r \ "-mm 

L 



£/r'^min (see 
According to 



(25) 
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Therefore, $(a;) should depend only on x and the ratio 
r/L. Indeed, Fig. ||b shows excellent data collapse for 
L/r — 8, with sharp cutoffs governed for a; < 1 by J\{x) 
and for x > {L/r)'^i^in by /2[a:(r/L)'^min]. 

In order to test the assumption that the functions /i 
and /2 are stretched exponentials with exponents 0^ and 
ipt, we plot 



nix)^log,o[A/^x)] 



(26) 



versus x in double logarithmic scale for various values of 
normalization constant A (see Fig. ^). If the stretched 
exponential conjecture is correct, Tl{x) should have two 
straight line asymptotes for logx — s- +oo with the slope 
tpg and for log a; — — cx) with the slope —(f>e. We find that 
the slopes <j>e and ipt^^^^^ of the straight line fits depend 
weakly on the value of A. Using A ~ 0.08, we obtain the 
longest regimes of straight line behavior. For this value 
of ^, we find « 2.1 and ipe « 2.5. Equation ( |l^ yields 
a predicted value of 0^ = 2.67 in good agreement with 
our simulation result. 



2. Off-critical behavior 

For p ^ pc, we identify three regimes determined by 
the value of the connectedness length, ^, in relation to 
the values of r and L: 

(i) ^ > L > 7-. In this regime, the fact that p ^ Pc can- 
not be detected because the connectedness length 
is larger than the other relevant variables. 

(ii) L> r. In this case, the upper cutoff of the dis- 
tribution Eq. (|l^) is governed by /a and the func- 
tional form of the rescaled probability <3> is given 

by 



(27) 



For large i, we suggest an exponential decay of <& 



$(£/r'*min) ~ exp ( -c 



(28) 



Indeed, for p < pc, semi-logarithmic plots of 
log $(£/r''miii) versus £ shown in Fig. |^a can be ap- 
proximated by straight lines with slopes which ap- 
proach zero as p ^ Pc- According to Eq. (p8|), these 
slopes k{p) should be proportional to ^^''min ^ 
\p — pf.\''-niin'^ « |p _ p^|i-i9. Figure IJb shows a 
double logarithmic plot of \k{p)\ versus \p — Pc\ for 
p < Pc- This curve can be well approximated by a 
straight line with slope 1.22 in good agreement with 
the scaling conjecture (|2^). For p > pc a similar 
analysis should hold. However, limitations on the 
size of the system we can simulate make the analy- 
sis problematic. Figure |[: shows P'(£) for various 



values of p > Pc- Note that it is only for values 
of p > Pc + 0.03 that the distributions "cut-off" at 
smaller £ than the distribution for p ^ pc- Thus it 
is only for values of p — pc > 0.03 that the large 
I behavior of Eq. ( p^ is determined by the fact 
that the system is not at criticality (i.e., by f^) as 
opposed to being determined by the finite size of 
the system (i.e., by /2). Below p = pc -I- 0.03, ^ 
is still greater than L. On the other hand, if p is 
not close to Pc, the scaling form is not expected to 
hold. Thus, the results are inconclusive based on 
the sizes of the systems we can generate — we can- 
not determine the parameters that govern the large 
i behavior of Eq. (ij) above Pc- 

(iii) L > r > When the connectedness length f 
is smaller than the distance r between the wells, 
the system can be considered homogeneous p8| , p7[ . 
This can be seen in Fig. ^ in which we plot P(i\r) 
for various values of r at p = 0.7 for 2d site per- 
colation (pc = 0.593). As r increases from below 
to above the connectedness length, the form of the 
distribution changes from the power law distribu- 
tion of Eq. ( p^ ) to a Gaussian distribution with a 
pronounced peak, a characteristic of homogeneous 
systems. Furthermore, as shown in Fig. pp, the 
fractal dimension of the shortest length crosses over 
from dmin = 1.13 to dmin = 1-0, characteristic of a 
homogeneous system [p8|j27|] . The convergence to 
a Gaussian can be expected due to the following 
considerations. The minimal path connecting the 
wells separated by distance r passes through r/^ in- 
dependent blobs. For each of these blobs, the prob- 
ability distribution for the shortest path across the 
blob, £b, is still given by Eq. ([T^), but with r and L 
replaced by ^ and £ replaced by £b. This distribu- 
tion is characterized by {£b) ~ ^'^min and variance 
0-2 ^ (^2^ _ (^j,)2 ^ ^2d„i„_ rpj^g minimal path 

is the sum of n — r/^ independent variables £b, 
hence it converges to a Gaussian with 



{£) - r^'^min-l 



and 



(j2 'p^'^^min ^ _ 



(29) 



Thus the slope of the graph, fc(p), of {£) vs r in 
Fig. 5c should decay as 



Hp) ^ |p-pcr''('*"""-i) = Ip-pcI""-'" (30) 



and the slope of cr^ versus r should decay as 



-''(2Vm-l) — 



(31) 



Indeed, (see Fig ||d) we see that the slope of cr^ 
versus r decays with p more strongly than that of 
{£) versus r. The numerical values of slopes from 
Figs. 1^ and |5|d are in good agreement with the 
theoretical predictions Eqs. ( |30| , |3l] ). For d = 3 we 
expect similar behavior. 
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D. Rectangular boundary conditions 



A. The model 



Since realistic oil fields do not have square boundaries, 
it is reasonable to ask what is the effect of rectangu- 
lar boundary conditions. Here we study the distribu- 
tions of the shortest length and minimal time on a two- 
dimensional lattice of size L^x Ly, ^ Ly. We position 
the wells at points A and B separated by distance r — 16 
along X-axis. We study two cases (Fig. ^): 

(a) Lx is fixed (L^ = 32), and we vary Ly {Ly = 64, 
128, 256, 512, 1024); 

(b) Ly is fixed {Ly = 32), and we vary L^ {L^ = 64, 
128, 256, 512, 1024). 

We find that (i) the shortest length and minimal time 
distributions are identical in all of the above cases and (ii) 
the scaling form for these anisotropic cases is the same 
as the isotropic case with L replaced by the minimum 
of Lx and Ly, i. e. the scaling form of the distribution 
Eq. (|8|) remains unchanged with the exception that L is 
determined by 

L = miii{Lx, Ly) . (32) 

Equation (|3^) can be a result of competing exponen- 
tials f2x and f2y. Since both f2x and f2y are rapidly 
decaying functions with L^ and Ly [see Eq. (1^)], the 
finite-size cutoff of P{£\r) is determined by the smaller 
of Lx and Ly. The fact that the results are independent 
of the axis along which the wells are aligned can be un- 
derstood by realizing that the probability that "oblong" 
clusters connect two points separated by a distance longer 
than the minimum of Lx and Ly is low. This finding has 
the implication that in anisotropic fields, a number of 
well pairs would be needed to optimize the recovery of 
oil in the field. 



III. MINIMAL TRAVELING TIME AND 
FASTEST PATH 

We turn next to dynamics, the study of flow on per- 
colation clusters, which has close ties to such applica- 
tions as h ydr ocarbon recovery and ground-water pollu- 
tion [Hj2^-^. In this section, we discuss the properties 
of the flow on d = 2 and d — 3 bond percolation clus- 
ters. Specifically, we investigate the scaling properties of 
the distributions of minimal traveling time and the length 
of the path corresponding to the minimal traveling time 
(fastest path) of the tracer particles. Some of the results 
in d = 2 were reported previously p6| . Here we extend 
the work to d = 3, and study the effects of finite system 
size and off-criticality for d ~ 2 and d = 3. 



We study incompressible flow between two sites A and 
B separated by Euclidean distance r (see Fig. 1). To 
model the flow front, we use passive tracers — particles 
not absorbed by the surroundings, that move only by 
convection, ignoring molecular diffusion which is slow on 
the time scales of interest. The convection is governed 
by the flow fleld due to the pressure difference between 
sites connected by the bonds. We simulate the flow of a 
tracer particle starting from the injection point A trav- 
eling through the medium along a path connected to the 
recovery point B. The dynamics of flow at a macroscopic 
level on the percolation cluster is determined by the lo- 
cal flow (local currents) on the individual bonds in the 
backbone of the cluster. The velocity of a tracer at each 
bond is determined by the pressure difference across that 
bond (Darcy's law ||33|]): 

V,, = T{P, - P.) , (33) 

where Pi and Pj are the values of pressure at sites i and j . 
The coefficient T, which is a function of permeability fc, 
viscosity rj, and the length of a bond Lb {T = k/{riLb)), is 
set to 1. We normalize the velocities assuming the total 
flow J between A and B is fixed, independent of the dis- 
tance between A and _B, and the realization of the porous 
media. This resembles more closely oil recovery processes 
where constant fiow, as opposed to constant pressure, is 
maintained. 

We obtain the pressure difference across each bond by 
solving Kirchhoff's law 

J2v.,^0, (34) 
j 

for each node i in the cluster where the summation is 
over all bonds connected to that node. Fig. |l]b shows 
the results of solving these equations for the cluster dis- 
cussed in Section I (Fig. |l|b). Magnitudes of currents 
on cluster backbone are depicted in gray scale with the 
lightest areas corresponding to the smallest currents and 
the darkest to the largest currents. 

We define the traveling time, t, of a path C as the sum 
of the tracer's traveling times at each bond {ij) joining 
sites i and j which are on the path, 

i= ^'1- (35) 
(«i)ec 

The traveling length, £, in turn, is the number of bonds 
present in path C. Among the ensemble of all paths {C}, 
we select the path C* that has the minimal traveling time, 

tmin: 

iinin(C*) mint(C) (36) 

{C} 

and we define the length of the fastest path, £niin, corre- 
sponding to the minimal traveling time, as the number 
of bonds present in path C*. The first quantity tmin is 
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the breakthrough time of the gas/hquid that displaces 
the oil during recovery and has fundamental importance 
to the oil industry. The quantity t determines post- 
breakthrough behavior. We also define the exponents 
dxj where x denotes ^mim iniin, ^ or t by 

X* ~ r'^" (37) 

and where x* is the characteristic (most probable) length 
or time of the corresponding distribution. 

Using a "burning" algorithm ||3^ , we then find the min- 
imal time and the fastest path for the particle to travel 
between points A and B. Figure [ij: shows the propaga- 
tion of the tracer particles through the same backbone 
shown in Fig. Qa and|^. At t = tmin, the tracer particles 
spread over t = <min • J bonds, which constitute a sub- 
set of the backbone with fractal dimension dtm, which is 
larger than the fractal dimension of the minimal path but 
smaller than the fractal dimension of the entire backbone 
de. Hence 

< rftm < ^B- (38) 



B. Minimal Traveling Time 

We first study the minimal traveling time for d = 2. 
In Fig. 0, a scatter plot of the minimal traveling time 
versus shortest path, we see that the minimal times are 
strongly correlated with the shortest paths in the real- 



izations simulated, tv, 



£^ , where z ~ 1.17. Since £ 



scales as r'^min we propose that tmin scales as r'^tm with 
dtm = zdjnin = 1-33. Tliis suggests that the same scaling 
form which applies to the distribution of shortest paths 
can also be applied to the distribution of minimal times, 
but with different exponents and amplitudes. Thus, we 
expect the Ansatz similar to Eq. ( p^ ) 



in 



fl 



h 



(39) 



where the scaling functions are /i (x) — exp(— atma^ ) , 
f2ix) = exp(-6tma;'^"°) and f^ix) = exp(-ctma;''tm). 
Here ^ is a characteristic length of the pair connected- 
ness function and has a power-law dependence on the 
occupancy probability p as 



C b - Pc 



(40) 



The first function /i accounts for the lower cut-off due 
to the constraint £ > r, while /2 and /a account for the 
upper cut-offs due to the finite size effect and due to the 
finite connectedness length, respectively. Either /2 and 
/a becomes irrelevant, depending on which of the two val- 
ues L or ^ is greater. For L < /2 dominates the upper 



cut-off, otherwise /s dominates. We have assumed inde- 
pendence of the finite size effect and off-criticality effect, 
so that Eq. ( |39| ) can be represented as a product of the 
terms which are responsible for the finite size effect (/2) 
and the effect of the concentration (/a). 

We sample over 10^ different realizations with the two 
sites A and B fixed. For each realization, we calculate ex- 
actly the minimal traveling time and the path which cor- 
responds to the minimal traveling time to obtain P(trnin) 
and P(Cin)- 



1. Behavior at criticality 

We first test numerically the scaling conjecture 
Eq. ( |39| ) at the percolation threshold p = Pc- In this case, 
^ = 00 and /a is a constant. Hence Eq. (p9) reduces to 



P'{tmin\r) 



1 /t ■ 

J- / ''mm \ 



{p^Pc). (41) 



Figure ||a shows that -P'(<mink) has a power-law regime 
with slope 



5^ = 2.0 ±0.1. 



(42) 



To determine the functions /i and /2, we compute the 
rescaled probability distribution 



P' (tmin |r) (imin)^^- r-'^tm (stm "D ^ (43) 



and plot it against scaling variable x = imin/''''*™ (see 
Fig. ||b). According to Eq. (|4^) 



'^{x)^Ah{x)h 



(44) 



Therefore, $(x) should depend only on x and the ratio 
r/L. Unlike the fractal dimension of the shortest path, 
dniin, there have been no calculations of the fractal dimen- 
sion of the minimal traveling time, dtm- We estimate dtm 
by finding the value which yields the best data collapse 
for Eq. (|^. For dtm = 1-33, Fig. ||b shows excellent 
data collapse with sharp cutoffs governed for small x <1 
by fi{x) and for large x > [L/rY^™ by f2[x{r/LY^^]. 

In order to test the assumption that the functions fi 
and /2 are stretched exponentials with exponents 0tm 
and i/'tm, we make a log-log plot of H(a;) = logjQ[^/$(a:)] 
versus x for various values of the normalization constant 
A (See Fig. ||c). If the stretched exponential conjecture 
is correct, H(a;) should have two straight line asymptotes 
for \ogx — > -f-oo with the slope i/'tm and for logo; —oo 
with the slope — (/)tm- The slopes 0tm and iptm of the 
straight line fits depend weakly on the value of A. Us- 
ing A — 0.14, we obtain the longest regimes of straight 
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line behavior. For this A we obtained (j)tm ~ 3.0 and 
"01111 ~ 3.0. With the same assumptions used to derive 
Eq. (n4), we can derive a similar expression for 0tm 



(45) 



1 



which yields a predicted value of (/)tm of 3.0 in good agree- 
ment with our simulation result. 



2. Off-Critical Behavior 

Finally, in order to test the dependence of -P'(tmink) 
on p we obtain data for large system size L {L = 1000) 
and for several values oip ^ Pc- As we do for the shortest 
length, we analyze the behavior of tmin in three regimes 
determined by the relation of the value of the connected- 
ness length, to the values of r and L. 

(i) ^ > L > r. In this regime, the fact that p Pc can- 
not be detected because the connectedness length 
is larger than the other relevant variables. 

(ii) L > ^ > r. In this case, the upper cutoff of the dis- 
tribution Eq. (|3^) is governed by /a and the func- 
tional form of the rescaled probability $ is given 
by 



(iii) L > r > ^. When the connectedness length is 
smaller than the distance between the wells, the 
behavior of the system is the same as a homoge- 
neous system [ P8y27| . This can be seen in Fig. [lO| a 
in which we plot P(imin|'') for various values of r at 
p — 0.6. As r increases from below the connected- 
ness length to above the connectedness length, the 
form of the distribution changes from the power law 
distribution of Eq. ( p| ) to a distribution with a pro- 
nounced peak, a characteristic of homogeneous sys- 
tems. In Fig. ^0|b, in order to eliminate the finite- 
size effect, we select L — r + 2 so that the distri- 
bution P(t\r) does not have a power-law regime, 
even for small r. In this case, as shown in Fig. |lO|c, 
the fractal dimension of the minimal traveling time 
crosses over from dtm = 1-33 to dtm — 2.0, charac- 
teristic of a homogeneous system ||2^,0 . The same 
considerations that we use to derive the behavior 
of the mean and variance of the shortest path can 
be applied to the mean and variance of the min- 
imal time. At the moment of breakthrough, i.e., 
when the first tracer particle reaches the second 
well, the part of the system filled with tracer parti- 
cles consists of rib = {r/O'^ independent blobs, each 
of which having a certain number of bonds {tmin)b 
with an average ((imin)&) = ^'^t™ and a variance 
Thus the average minimal time for the 



entire system scales as 



h 



(46) 



For large t„iin, we suggest an exponential decay of 



$(imin/r''t'°) - exp 



in 



(47) 



Semi-logarithmic plots of log <i>(iinin/?''^'™) versus 
trnin for p > pc and p < Pc shown in Fig. ^ and 
yb, respectively, can be approximated by straight 
lines with slopes which approach zero as p — > Pc- 
According to Eq. (|47|), this slope k{p) should be 
proportional to 



(^min) — 



rib 



with a variance 



(49) 



(50) 



The scahng plot (Fig. [l0|d) of (tmin) versus \p — pc\ 
shows good agreement with the theoretical predic- 
tion of Eq. (11 



= iP-Pcr' 1 



d = 2]. 



(51) 



\P-Pc 



\P~Pc 



11.77 



(48) 



Figure ||c shows double logarithmic plots of \k{p)\ 
versus \p — Pc\ for p < Pc and p > Pc, which can 
be well approximated by straight lines with slopes 
1.81 and 1.77, respectively, in good agreement with 
the scaling conjecture, Eq. (p8|). As was the case 
with the analysis of P'{£\r) > Pc iov d = 3 (see 
Sec. II.C.2.ii), we cannot determine the parameters 
that govern the large tmin behavior of P'(imin) be- 
cause of limitations on the size of the system we 
can simulate. 



The graph of a versus r (see Fig. |10|e) shows linear 
behavior of a versus r in agreement with Eq. ([S^). 
Equation (^ also predicts that the slope of this 
linear dependence decays as 



\p-Pc 



-[dtm~(d/2)]i^ 



\P-Pc 



-0.42 



[d = 2]. 
(52) 



However the measured slope has very small varia- 
tion with \p — Pc\ which is beyond the accuracy of 
our data points. 
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As mentioned above, the minimal traveling time is the 
sum of the inverse local velocities over the fastest path 
where the fastest path is statistically identical to the 
shortest path. While the velocity distribution has been 
studied extensively |^,^, because the velocities along 
the path are correlated, the relation between the min- 
imum traveling time distribution and the local velocity 
distribution is an open challenge for further research. 

The analysis for three dimensions is completely analo- 
gous to that for two dimensions. Our results are shown 
in Figs. |ll| and |l^ and the scaling parameters found are 
included in Table I. 



C. Fastest Path 

We observe that the path which takes minimal time is 
not always the shortest path. However analysis of the dis- 
tributions of ^„iin yields parameters identical to those for 
the distribution of the shortest paths between poi nts sep- 
arated by distance r studied in detail in Ref. [|5|. Thus, 
statistically, the path which takes the shortest time is one 
of the paths of shortest length. 

In many transport problems, the characteristic time t* 
scales with the characteristic length £* with a power law, 

t* - {tf . (53) 

Since t* scales as r'^' and £* scales as r'^min , 

z = dt/d^in. (54) 

Since tmin Q-nci ^min 

are strongly correlated, the distribu- 
tions P(^min) and P(tmin) Satisfy 

P(^min)^^inin — P(^min)^^inin (55) 

Combining Eqs. (^3|)-(^) and the equations for the re- 
spective distributions, we obtain a scaling relation be- 
tween exponents, 

(5Vin ~ l)^Vin (ftm - l)dtm- (56) 

These scaling relations are well satisfied by the set of 
scaling exponents given in Table I. 

IV. CONCLUSIONS 

By modeling porous media by bond percolation and 
using concepts of percolation theory, we study the flow 
of fluid in porous media in 2 and 3 dimensions between 
two "wells" separated by Euclidean distance r. We in- 
vestigate the distribution function of the shortest path 
connecting the two sites, and propose a scaling Ansatz 
that accounts for the dependence of this distribution (i) 
on L, the size of the system, and (ii) on p, the bond occu- 
pancy probability. The finite size of the system, L corre- 
sponds to the size of the oil field and the bond occupancy 



probability, p is related to the good sand to impermeable 
rock ratio of the porous media. We confirm by extensive 
simulations that the Ansatz holds for d = 2,3, and we 
calculate the relevant scaling parameters. 

In order to understand the properties of the flow of oil 
displaced by fluid or gas, we study the dynamics of flow 
on percolation clusters. We study two dynamical quan- 
tities: the minimal traveling time and the length of the 
path corresponding to the minimal traveling time. Be- 
cause of the approximate parallel between the shortest 
path and the minimal traveling time of flow, the study 
of the shortest path is the first step in understanding the 
properties of the oil fields. In particular, a scaling Ansatz 
for these dynamical quantities includes the effect of finite 
system size and off-critical bond occupation probability. 
We find that the scaling form for the distribution func- 
tions for these dynamical quantities for d = 2, 3 is similar 
to, but not identical to, that for the shortest path. In ad- 
dition to calculating the relevant distribution functions 
and scaling relations, we also calculate a number of new 
exponents (see Table I). 
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APPENDIX A: SAMPLING ALL SIZE CLUSTERS 
VERSUS SAMPLING ONLY INFINITE 
CLUSTERS 

Since P'{£\r) is the distribution of shortest paths be- 
tween point in clusters of all sizes (not just the infinite 
cluster) , we can write 

/•OO 

P'{£\r)^ I P{£\r,M)P{M)dM, (Al) 
Jo 

where P{£\r, M) is the distribution of shortest paths in 
clusters of mass M. 

The size L of a cluster of mass M is of the order 
jV^i/d/ and the shortest path between two points in the 
cluster is larger than L'^mhi drops off exponentially, so 
P{£\r, M) will effectively be zero when M is less than 
^'^//'^min. On the other hand, when M is greater than 
this value, P{l\r, M) will be the same distribution as the 
distribution for the case when the two points are in an in- 
finite cluster. This can be taken into account by replacing 
P{£\r, M) with P{£\r) and replacing the zero lower limit 
of the integral by £'^//Vin. Since P{M) ~ M""', where 
T = d/df, we have 




P(£\r)M-'^''^fdM 

'mill 



= P(£|r)^('*/-'*)/''nnn^ (A2) 
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from which follows 

9't = 9i + ^^^=9i + {d-df)i). (A3) 

This can be generalized for any quantity, x, with fractal 
dimension dx and with distributions 

nAr) = ^^(;^y" (A4) 

and 

where "infinite clusters" and all finite-sized clusters, re- 
spectively, are sampled. This generalization results in 

5x = 5x + ^^- (A6) 
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Table I. Summary of exponents and coefficients in scaling form 

P{x\r) ~ {l/r'^-){x/r'^-)-^''Mx/r'^-)f2{x/L^-)fs{x/^'^-) 

where fi{y) = ex.p{—axy~'^''), fiiv) = e^pi—bxy^"), fsiv) = exp(— c^y). Here x denotes one of the quantities, i or 
t,nin- The notation N/A means Not Applicable (since no theoretical value exists), while the notation (+/-) indicates 
above or below Pc- 

(a) d = 2 results 



X 

exponent 


I 




Sim. 


th(X)rv 


Sim. 


th{X)ry 




1.13 ± 0.01 


.V/.l 


1.33 ± 0.05 


N/A 


g'x 


2.14 ±0.02 


2.11 


2.0 ±0.1 


N/A 


0;,. 


0.5 


NjA 


1.1 


N/A 


(:>j. 


7. .3 ± 0.5 


l/((/,, - 1) = 7.()9 


3.0 


3.0 


bx 


3.5 


NjA 


5.0 


N/A 




4.0 ±0.5 


N/A 


3.0 


N/A 


C-x 


2.4(-),3.7(+) 


N/A 


1.6(-),2.6(+) 


N/A 



(b) d = 3 results 



X 

exponent 


i 


^min 


sim. 


theory 


sim. 


theory 


dx 


1.39 ±0.05 


N/A 


1.45 ±0.1 


N/A 


9'x 


2.3 ±0.1 


2.23 


2.1 ±0.1 


N/A 


ax 


1.4 


N/A 


2.5 


N/A 


4>x 


2.1 ±0.5 


l/(d, - 1) - 2.56 


1.6 


2.0 


bx 


2.0 


N/A 


2.3 


N/A 


4'x 


2.5 ±0.5 


N/A 


2.0 


N/A 


C-x 


3.1(-) 


N/A 


2.9(-) 


N/A 
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FIG. 1. An example of (a) percolation backbone (thin lines) and shortest path (thick line) corresponding to the flow between 
points A and B — note that the lines here represent not microscopic pores but rather sand bodies whose size are of the order of 
tens of meters; (b) magnitudes of currents on cluster backbone are depicted in gray scale with the lightest areas corresponding 
to the smallest currents and the darkest to the largest currents; and (c) time evolution of the flow between points A and B 
on a two-dimensional percolation cluster. Thick lines denote bonds reached by tracer particles at time tmin- The double line 
between A and B denotes the "fastest path" between these two points. 
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FIG. 2. (a) Log-log plot of P(r\£) for two-dimensional percolation at the percolation threshold for ^ = 64 (O); 256 (□) 
and 1024 (A). The best data collapse is obtained with = 0.88 and the slope of the power regime is 2.11. (b) Log-log plot 
of P{l\r) (+) and P'{l\r) (Q) for two-dimensional percolation at criticality and for the system size L = 1024 and the distance 
between wells is r = 32. The power-law regime of P' {l\r) has slope g'l = 2.02 (thick solid line), while that of P{t\r) has slope 
ge = 1.95 (thin solid line). The purpose of this figure is to illustrate the difference between gt and g'^. The values shown here 
are lower than those predicted by Ziff |^ because we have not included corrections to scaling in our determination of these 
quantities. 
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FIG. 3. For d = 3, (a) log- log plot of P'{£\r) at criticality {p = Pc ~ 0.2488) and for different sets of parameters: 
{r,L) = (2, 32), (4, 64), (8, 128). The straight line regime has slope g'^ — 2.3. (b) Log-log plot of rescaled probability 
${x) = P'{£\r) x^t r'^min against rescaled length x = (./r'^min using the values, g'^ — 2.3 and dmin = 1-39. The curves 
are flat in the center because f2{x) is a stretched exponential (see Eq. (§^). (c) Log- log plot of transformed probability 
n(a;) = login [A/$(a:)1 versus x — £/r''min. The slopes of the solid lines give the power of the stretched exponential function 
/i and /2 in Eq. (|25[). Using the parameter A = 0.08, the slopes give ~ 2.1 for the lower cut-off and ip ~ 2.5 for the upper 
cut-off. 
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FIG. 4. For d — 3, (a) semi-logarithmic plot of transformed probability ^(l) (see Eq. (pq)) versus £ shows pure exponential 
behavior of /a. (b) The slope of the log-log plot of the coefficient in exponential function as a function of |p — Pc| gives the 
value i^dmin ~ 1-22 for p < pc- (c) P'{£) for p > pc- Note that it is only for p > 0.2788 that the large £ behavior is determined 
by the fact that the system is not at criticality. 
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FIG. 5. (a) Distributions of P{£\r) for r = 4, 8, 16, 32, 64, 128, 256, 512 and for p = 0.7. To reduce the lattice effects, data 
is obtained for the pairs of wells on the a;-axis. (Note that for this case, where r > ^, the distributions P'{£\r) and P{£\r) are 
essentially the same since all the clusters span the lattice.) The distributions converge for large r to a Gaussian with mean {£) 
shown in parts (b,c) and variance = {^^) — {£)^ , shown in part (d) as functions of r for p = 0.65(0) P ~ 0.7(a). (b) 
Log-log plot of {£) versus r. Note the crossover from power law behavior with exponent d,niii ~ 1-13 to linear behavior with 
exponent 1.0. (c) Same as (b) on linear scale. The slopes of linear the fits k{p) are 1.45 for p — 0.65 and 1.30 for p = 0.7. This 
yields k{p) ~ (p — Pc)~°'^^ in good agreement with equation Eg. (|3^) . (d) The dependence of versus r. According to Eq.(p9[), 
the dependence becomes linear only for r > ^ ~ (p — pc)~" , indicated on the graph. The slopes of linear fits k{p) are 0.33 for 
p = 0.65 and 0.12 for p = 0.7. This gives k{p) ~ (p — Pc)^^'^ in good agreement with equation Eq.(pl[). 
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FIG. 6. (a) Two configurations with rectangular boundaries are studied. (b) For r = 4, different system sizes 
{L^,Ly) = (32, 128), (32, 256), (32, 512), (32, 1024), (128, 32), (256, 32), (512, 32) and (1024,32) are studied and all the distri- 
butions of minimal path are collapsed into each other. 
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FIG. 7. For d = 2, scatter plot of the minimal traveling time tmin versus minimal length £ and for a fixed well separation 
r = 1. Note the strong correlation between tjnin and i. The slope of the tail of the scatter plot is 1.17 yielding a value of 
dim. = 1.17 • dmin = 1.32, consistent with our result in Table I. 
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FIG. 8. For d = 2, (a) log-log plot of P'{t\r) for p = Pc = 0.5 and for different sets of parameters, 
(r, L) = (16, 250), (32, 500), (64, 1000). The straight line regime has slope g't = 2.0. (b) Log-log plot of rescaled probability 
${x) = P' {tj^-iiji\r)x^tr'^* against rescaled length x = tmip /r'^' using the values, g't = 2.0 and dt = 1.33. The curves are flat in the 
center because f2{x) is stretched exponential (see Eq. (P5|)). (c) Log-log plot of transformed probability n(a;) — loRTn\A/^(x)] 
versus x = ^min/^'**- The slopes of the solid lines give the power of the stretched exponential function /i and /2 in Eq. (E5|). 
Using the parameter A = 0.14, the slopes give (j> « 3.0 for the lower cut-off and ~ 3.0 for the upper cut-off. 
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FIG. 9. For d — 2, (a) somi-logarithmic plot of transformed probability $(tmin/^''') versus tmin for /s for 
p = 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48 below criticality. (b) Same for p = 0.52, 0.53, 0.54, 0.55, 0.56 above criticality. (c) 
The slope of the log-log plot of the coefficient in exponential function /a as a function of |p — Pc| gives the value udt w 1.77 for 
p > Pc and 1.81 for p < pc. 
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FIG. 10. (a) Log-log plot of P(tmin |r ) for p = 0.6 and for r = 4, 8, 16, 32, 64, 128, 256 and L = 258. The distributions for large 
r converge to Gaussians with mean (tmin) and variance a^. (b) Log-log plot of P(tniink) for p = 0.6, r = 4, 8, 16, 32, 64, 128, 256 
and L = r + 2. (Note that for this case, where r > 5, the distributions P'(tmin|'') and J-'(tminl^) are essentially the same since 
all the clusters span the lattice.) (c) Log-log plots of (tmin) versus r for p = 0.6 and L = r + 2. (d) Log-log plot of the scaled 
average minimal traveling time, (tmin)/'''^, versus p — Pc for r = 128, 192, 256, 384, 512 and L = r + 2. Note that in all cases 
r ^ The slope of the line, 0.84, is in good agreement with the theoretical prediction, 0.89. (e) The behavior of the width, a, 
of the distributions of the travcliug time versus r for p = 0.53, 0.54, 0.55, 0.57 and 0.6. The graph shows approximately linear 
dependence of a on r. The variation of the slope with p — Pc is within the error bars of the data. 
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FIG. 11. For d = 3, (a) log-log plot of P'{t\r) for p = pc = 0.2488 and for different sets of parameters, 
(r, L) — (4, 32), (8, 64), (16, 128). The power-law regime has slope g't — 2.1. (b) Log-log plot of rescaled probability 
"I>(2;) = P'(tniin|'")a^^'^'''' against rescaled length x — tmip /r'^' using the values, g'^ — 2.1 and dt — 1.45. The curves are flat in the 
center because f2{x) is stretched exponential (see Eq. (P5|)). (c) Log-log plot of transformed probability n(a:) — loginf^/^(a^)1 
versus x = imin/^'*'- The slopes of the solid lines give the power of the stretched exponential function /i and /2 in Eq. (E5|). 
Using the parameter A — 0.08, the slopes give (j> ~ 1.6 for the lower cut-off and ip ~ 2.0 for the upper cut-off. 
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FIG. 12. For d = 3, (a) semi-logarithmic plot of transformed probability ^'(tinin) versus tmin below critical point for 
p = 0.1988, 0.2088, 0.2188, 0.2288.0.2358, 0.2388 shows pure exponential behavior of /s. (b) The slope of the log-log plot of the 

coefficient in exponential function fs as a function of \p — pc\ gives the value i^dt « 1.30 for p < pc- (c) P(tmin) for p > Pc- 
Note that for the values of p simulated, the large tmin behavior is determined by the finite size of the system — not /s. 
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